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Abstract 

We have exploited a variety of techniques to study the universaHty and 
stabihty of the scahng properties of Harper's equation, the equation for a 
particle moving on a tight-binding square lattice in the presence of a gauge 
field, when coupling to next nearest sites is added. We find, from numerical 
and analytical studies, that the scaling behavior of the total width of the 
spectrum and the multifractal nature of the spectrum are unchanged, provided 
the next nearest neighbor coupling terms are below a certain threshold value. 
The full square symmetry of the Hamiltonian is not required for criticality, 
but the square diagonals should remain as reflection lines. A bicritical line 
is found at the boundary between the region in which the nearest neighbor 
terms dominate and the region in which the next nearest neighbor terms 



dominate. On the bicritical line a different critical exponent for the width of 

the spectrum and different multifractal behavior are found. In the region in 

which the next nearest neighbor terms dominate the behavior is still critical 

if the Hamiltonian is invariant under reflection in the directions parallel to 

the sides of the square, but a new length scale enters, and the behavior is no 

longer universal but shows strongly oscillatory behavior. For a flux per unit 

cell equal to 1/q the measure of the spectrum is proportional to 1/q in this 

region, but if it is a ratio of Fibonacci numbers the measure decreases with a 

rather higher inverse power of the denominator. 
73.40.Dx,73.50.-h,03.65.Sq 
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I. INTRODUCTION 



Harper's equation can be derived as the equation for an electron in a strong two- 
dimensional periodic potential and a weak magnetic field, or for an electron in a strong 
magnetic field and a weak periodic potential. The Hamiltonian can be written in the form 

H{px,x) = 2ta cos px + 2tb cos X , (1.1) 

where 

Px = —2ni(pd/dx (1.2) 

is the variable conjugate to x. AzbelS showed that in the case ta = tb, which corresponds 
to a periodic potential with square symmetry, the spectrum forms a "devil's staircase" for 
irrational values of 0, and Hofstadteii generated computer drawings of the spectrum for 
rational values of 0, and discussed the self-similarity and scaling of the spectrum. Work by 
Aubry and Andre@ exploited the symmetry of the problem under canonical transformations 
that interchange x and px, and showed that, for irrational values of 0, there is a localization 
length in the x direction, independent of energy, 

L = l/ln(V4), for tb > ta, (1.3) 

which diverges at this symmetry point ta = tb- Also the sum of the widths of all the energy 
bands, the measure of the spectrum, has the form 

W = A\tb-ta\, (1.4) 

which vanishes at the same point. These properties have suggested that this point is like a 
critical point of the system, that |1 — ta/tb\ represents the distance from the critical point, 
and that, where = p/g is a rational, the denominator q acts as a finite size in finite size 
scaling theoryi. 

Various methods have been used to study the Harper equation. There are some rigorous 
analytical results on the Lyapunov exponent (reciprocal of the localization length)!, and on 



the measure of the spectrum (sum of the band widths)Qilll. Aubry duahty gives information 
about the locahzation length. 

If H{px,x) defined in Eq. ( |1 . Ij ) is treated as a classical Hamiltonian there is an obvious 
interpretation of the critical point, since for ta > tb the energy contours surrounding the 
minima of H are separated from the contours surrounding the maxima of by a region of 
orbits open in the x direction, while for ta < tb there are orbits open in the px direction. For 
the symmetric case ta = tb there is a single energy E = separating the two types of closed 
orbits. A more subtle semiclassical analysis, based on scaling theory, must be introduced to 
explain why the behavior is not just singular at the energy of this separatrix, but is singular 
at all energies in the spectrum. 

Much of the information we have about this problem comes from numerical studies of the 
spectrum for rational values of similar to the one carried out by Hofstadter. The vanishing 
of the measure of the spectrum at the critical point gives a one-parameter indicator of the 
critical point. Numerical studies on the Harper equation indicate that this measure is given 
by Eq. ( |1.4| ) for all irrational values of (f), and that for rational values (p = p/q the measure 



where the scaling function g{s) behaves as 9.3299/s when its argument is small, and tends 
to ±4 for large Only the corrections to scaling seem to be sensitive to the value of jJUll. 
The distribution of band widths within this spectrum gives a more detailed picture of the 
spectrum, but this is sensitive to the value of 0. The simplest forms might be expected if 
is an irrational solution of a quadratic equation, so that its continued fraction expansion 
repeats itself, or if it is a rational approximant obtained by truncating such a repeating 
continued fraction. The golden mean, which is the limit of the ratios of successive Fibonacci 
numbers, is the simplest case of this sort. In this case the spectrum has a multifractal 
structure which is self-similar at the critical pointiiHli. 




Another method which has been used to get a number of important results takes its 



scales ai 



W-{tb- ta)g[q{l - ta/tb)] , 
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simplest form when is small, so that Eq. ( p..l| ) can be treated semiclassically. It is then 
a second order difference equation with a slowly varying central term, for which the WKB 
equation can be used. For ^ ta only those bands close to E = have appreciable band 
width, and the others all have widths which vanish exponentially as 1/0 gets large, since 
the widths are determined by tunnelling between successive minima or maxima of 2t;,cosa:. 
The case = 1/g is particularly simple, and this has been used to derive an analytic form 
for the function g{s) of Eq. (pT5|)i'ii. 

The critical properties of the Harper equation have been explored in great detail by 
various combinations of these methods. However, it is not clear to what extent the elegant 
properties of the Harper equation are special properties of that equation, and to what extent 
they are robust, and survive perturbations of the model. On the basis of the small behav- 
ior, Suslov§ has argued that modifications of the x-dependence of Eq. ( |1.1|) which maintain 
the periodicity will lead to an energy-dependence of the critical value of tt- This behavior is 
supported by numerical calculations.EUll On the other hand Helffer and Sjostrand0 have 
argued that variants of Eq. ( |1 . 1| ) which preserve the invariance under the canonical trans- 
formation critical at all energies. It is the purpose of this work to take 
these questions further by exploring in some detail a simple generalization of the Harper 
equation. 

The equation we explore is the generalization of Eq. (|1 . 1| ) which takes the form 

H{px, x) = 2ta COS + 2tf, COS X + 2t„^ cos{px — x) + 2tab cos{px + x) , (1.6) 

This could represent a tight-binding model in which electrons can tunnel to next nearest 
neighbors as well as to nearest neighbors on a rectangular lattice, and, for tab = 0, = 
t{, = t^i, it can represent an electron on a triangular lattice^. This form of the Hamiltonian 
is very convenient for numerical work, as it remains in the form of a three-term difference 
equation, and there have been a number of earlier studies of itSSIl. It has enough parameters 
that we can study the effects of different features of the Hamiltonian such as the breaking 
of the symmetry, and the change in the nature of the constant energy contours. For ta = tb 
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and tab = tab the system has square symmetry, but there are two quite different regimes with 
this symmetry, according to whether ta or tab is dominant. We show that there is actually 
an interesting bicritical point which separates the two different regimes. 

In this paper we have exploited most of the techniques mentioned above. In Sec. II we 
discuss the main features of the energy contours of Eq. and discuss the symmetries 

of the system. In Sec. Ill we discuss the form of the characteristic equation whose roots 
give the eigenvalues, repeat some earlier results which were obtained by the use of Aubry 
duality, discuss the extension of finite size scaling arguments to this case, and extend the 
sum rule of Last and Wilkinson0 to this case. In Sec. IV we present an analysis in terms 
of multifractals for the case that is a ratio of neighboring Fibonacci numbers. In Sec. 
V we give the results of calculations of the measure of the spectrum when is a fraction 
such as 1/q or p/{p'^ + 1) that represents a slow modulation of the diagonal term of the 
difference equation. These two approaches complement one another, since for the ratio of 
Fibonacci numbers the spectrum is spread over a large number of bands, even close to the 
critical point, whereas for small values of (f) the only bands with significant measure are 
close to zero energy. We find that the behavior is relatively simple in the region in which 
ta = tb is dominant, but there is a bicritical region not only at ta = tb = 2tab = '^tab, but 
also for ta = 2tab = ^tab > tb- For tai = tab dominant the situation appears to be much more 
complicated, with some important oscillatory terms which confuse the analysis of numerical 
results. 

In Sec. VI we do what we can to explain the results we have from numerical analysis 
in terms of WKB theory and scaling theory. In some cases our understanding is reason- 
ably complete, but in other cases we can do little more than explain why the problem is 
complicated. There is a concluding discussion in Sec. VII. 
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II. CLASSICAL ORBITS AND SYMMETRY 



In this discussion we assume all the coefficients ta, h, t^j,, tab are positive. For the case 
ta = tb, t^i = tab the spectrum of the Hamiltonian given by Eq. ( |1.6| ) is invariant under 
the eight operations of the symmetry group of the square. The four proper rotations are 
generated by —>■ —x, x ^ Px, while the time reversal operation p^ —p^ generates the 
improper rotations. The Hamiltonian is also invariant under the group of translations in 
phase space corresponding to a square lattice. The classical Hamiltonian has a maximum 
at Px = = X, where its value is 4(ta + tab), and at equivalent lattice points. For ta > 2tab it 
has a minimum at Px = n = x, where its value is —4:(ta — tab)- There are two saddle points 
in each unit cell, at px = 0, x = vr, and at px = 7^, x = 0, where its value is —4:tab- At this 
value of the energy there is a contour given by 

2tab + ta cos X 

COSPx = (2.1) 

ta + 2tab COS X 

which threads the system, separating contours that surround minima from those that sur- 
round maxima. For ta < 2tab the maximum is unchanged, but the points at Px = it = x 
become subsidiary maxima, while the points at px = 0, x = vr, and at px = tt, x = 
become minima. Four more saddle points appear at cospx = —ta/2tab, cosx = —ta/2tab, 
where the energy is —t^/tab- The contour joining the saddle points and separating contours 
surrounding minima from those surrounding maxima is now 

cos Px = -ta/2tab OT COS X = -ta/2tab ■ (2.2) 

In this paper we pay particular attention to the bicritical point ta = tb = 2tab = 2tab, where 
all contours surround maxima except for the lines where the energy has its minimum value 
—2ta. At the points px = it = x the lowest nonvanishing partial derivatives of the energy 
are the fourth derivatives. 

For the case ta = tb, tai > tab the symmetry operation px x, a reflection symmetry 
in phase space, remains, as well as rotation by vr. There is still a maximum at px = = x, 
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where the energy is ita + 2tab + '^tab, and at equivalent lattice points. For ta > tai there is 
a minimum at = = x, where its value is — 4ta + 2t^5 + 2tab- For > 4:t^itab there are 
two saddle points in each unit cell, at Px = 0, x = tt and at Px = tc, x = 0, where its value 
is — 2(t„5 + tab)- At this value of the energy there is a contour given by 

-1 



2tabCOs[-{px + X)] = -COs[-{px - X)][ta± ^Jtl- At^ltab] (2.3) 

which threads the system, provided ta > t^i + tab- For ta < tab + tab there is a range of 
energies in the neighborhood of —2{tab + tab), with cos[^{px — x)] near ±1, for which there 
are no values of [px + x) that satisfy Eq. (|2.3| ). In this case there are open orbits in the 
direction of constant {px — x). 

A special case of this sort is ta = tb = tai, tab = 0, which has triangular symmetry, and 
is equivalent to the case worked out numerically by Claro and Wannierill. 

For ta tb, t^i = tab the Hamiltonian is invariant under px — > —px- For tab = tab > tb > ta 
there are maxima aX px = = x and a.t px = it = x, where the energy has the values 
4:tab ± (ta + tb), and minima at px = 0, x = n and at Px = tt, x = 0, where the energies are 
—4:tab ± (ta — th). There are saddle points given by 

cos = -tb/2tab, COSX = -ta/2tab , (2.4) 

and the energy has the value —tatb/tab on the lines on which either of these conditions is 
satisfied. For tb > 2t^i = 2tab > ta the four saddle points given by Eq. ( |2.4| ) disappear, and 
the new saddle points are at px = 'n', x = and at px = tt = x, where the energies are 
—ta ± {2tab — tb). There are orbits open in the px direction between these energies. 

III. CHARACTERISTIC EQUATION AND AUBRY DUALITY 



From Eq. ( |1.2| ) it can be seen that the operator 2cosp^ is a displacement operator that 
displaces the coordinate by 2770. The eigenvalue problem for the Hamiltonian ( |1.6| ) takes 
the form of a set of finite difference problems 

(ta + t^-^e2-^(""5)+*fc2 + ^^^^-2n^Hn-^^)-^k,^^^_^ ^ 2^ COs(27r0n + ^2)0^ 



+ {ta + t^5e-2-'^(-+5)-^^2 + Ue2-'^("+^)+^^^)a„+i = Ea^ , (3.1) 

with the variable parameter ^2 determined by the initial value of x. When = p/g is 
rational this equation is periodic with period q, and solutions of the Floquet form 

a„ = c^e"'^^ , (3.2) 

with c„ periodic, can be found. This then yields a finite matrix problem, with the matrix 
tridiagonal apart from the top right and bottom left corners, for which the characteristic 
polynomial has as its only ki dependent term 

q-l 

^_iy-i^ik^g Yl (ta + Vbe"^^"^^^"'''' + Ue^("+^)+'^^) + conjugate complex . (3.3) 

ra=0 

This product can be expressed in terms of a Chebyshev polynomial Tg (see Appendix A), 
and this gives 

+ (-l)^2{t^5 cos[g(A;i - fcs)] + cos[q{h + ^2)]} , (3.4) 

where 

= iab^ab ■ (3.5) 

Since the whole spectrum is, from Eq. (|1.6| ), clearly invariant under the interchange of p^, x 
and ta, tb, there must also be a similar ^2 dependent term in the characteristic polynomial, 
so the characteristic polynomial can be written in the form 

P{E) = Po{E) - {-iy4t^{cos{qh)T/^^) + cos{qk,)T/±)} 

+ (-l)P2{4cos[g(A:i - k,)] + tl,cos[q{ki + k2)]} , (3.6) 

where Po{E) is independent of ki, k2- The energy bands are determined by the variation of 
the solutions of the characteristic equation as ki, k2 are varied. 
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Some simple analysis can show which of the terms in this expression will dominate in 
the limit of large q. For a; > 1 we have 

[Tq{x)Y''^ ^ X + , (3.7) 

and so, for tb > ta, t^i > tab, the dominant term in Eq. ( p. 61 ) is of order 

[litf, + ^/tl-4t^)Y (3.8) 

for tb > t^i + tah, and is of order for tb < t^i + tab- This transition from a regime where the 
band widths are dominated by the term depending on k2 to a regime where the band widths 
are dominated by a term depending on ki — k2 occurs at the same values of the parameters 
as the changes in the nature of the classical orbits which we discussed in connection with 
Eqs. ( p.3|) and (|2.4| ). Under these conditions only one term in Eq. ( p.6|) is relevant for large 
q, or two terms if ta = tb- Since this expression is independent of the value of the energy, one 
should expect the critical values of the parameters to be energy independent, whereas the 
classical orbits only give information about the behavior at the singular value of the energy. 

For 2tab = '^tab > tb>ta the situation is very different, since the Chebyshev polynomials 
are now of order unity, and all four terms in Eq. ( p. 61 ) would seem to be marginal for large 
q. In particular one should expect the band widths to depend on 

cos[g arccos(ta/2ta5)] and cos[qaxccos{tb/2tab)\ , (3.9) 

so the band widths should display nearly periodic behavior in q. In fact we found such 
behavior in the numerical studies reported in this paper before we had realised that they 
ought to be found. 

In an earlier worki it was shown how the argument of Aubry and Andre can be adapted 
to this situation. There are three parts to this argument. Firstly they state that the element 
of the Green function connecting the two ends of a tridiagonal matrix can be expressed as the 
product of the next-to-diagonal matrix elements divided by the characteristic polynomial. 
The product of off-diagonal matrix elements is just the coefficient of cospx in Eq. ( p. 61) . 
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Secondly they compare this expression with the expression for the dual problem obtained 
by interchanging and x. The characteristic polynomial is unchanged, and the ratio of 
the products of next-to-diagonal elements can be used to generate an expression for the 
difference between the Lyapunov exponents in the x and px directions, in the form 

This is zero for t^i + tab >tb>ta. For tb>ta> t^i + tab it gives 



2 -a I 2 

and for > 1^1 + tab > 4, tab > ^ab it gives 



A. - A, = ln[i^^±MLif!] , (3.11) 



A. - A, = ln[ ^ 7 " ] . (3.12) 

tab 

The third part of the argument, which we find to be rather more subtle, says that if A^, 
is positive then Ap must be zero, so Eqs. ( |3.11|) and ( p.l2|) are actually equations for A^, 



rather than for A^; — Ap. Eigenstates for the Hamiltonian (|1.6| ) in px space can be found 
from eigenstates in x space by Fourier transformation. These eigenstates in x space have 
their support on a lattice of points, so their Fourier transforms are periodic. If they are 
localized in space their Fourier transforms are smooth periodic functions. Functions of this 
sort corresponding to the same value of the energy cannot be superposed to give localization 
in Px space as well as localization in x space. 

These results show that states are localized in x space, independent of energy, for generic 
irrational provided tb is greater than both ta and tai + tab- This condition, now independent 
of energy, is the same as the condition for the existence of open orbits extended in the px 
direction given in the discussion of Eq. (|2.4|) . 

The finite size scaling argument that lead to Eq. ( |1.5|) can be generalized to deal with 
the critical properties of Eq. ( |1.6|) . For tb > ta > tab + tab the width for irrational (p is still 



given! by Eq. ( |1.4| ) , but the scaling length is now given by Eq. ( p.ll|) so that we have 



\tb + k\tl-At^ tb-t. 



w>,itb- ta)g{qH l' I /I .J ) ^ - iMq^===) ■ (3.13) 

¥a + ytl-At^ Jtl - 4t2 
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In particular, this results in the prediction that the measure of the spectrum should scale as 

9.3299 



X — -— ^t2_4^2 (314) 

for ta = h > t^i + tab- A special case of this is the triangular lattice with ta = h = t^ii 
tab = 0, where the measure of the spectrum is the same as it is for the square lattice. 

For the case tf, > 2tai = ^tab > ta the measure of the spectrum for irrational is itb — Stab 
and the scaling length is given by Eq. (|3.12|) , so the use of the same argument would lead 



to a finite size scaling expression of the form 



-\/^b + 2tab + v^tfc — 2t, 



W X (t, - 2t)G{qVh^b '' ' "^"'' ) . (3.15) 

^tab 



The scaling function G{s) must diverge as s at the origin, so that at the point tb = 2t 



ab 



2tab the measure of the spectrum is finite, and goes to zero like for large q. However, this 
argument does not take account of the second length scale introduced in Eq. ( |3.9| ), so we 
should expect the function G in Eq. ( |3.15| ), and the coefficient of the limiting behavior. 



to depend on ta according to the form given in Eq. (|3.9| ), which is periodic or nearly periodic 
in q. 

The argument of Last and Wilkinson^ which provides a lower bound to the spectrum 
for the critical case can be generalized to deal with the situations we consider in this work. 
The simplest case is given by tb > ta > tai + tab, where the result of Avron, Mouche and 
Simoni that the intersection spectrum (the intersection over ^2 of the spectra for fixed 
has measure 4(tfe — ta) remains valid, as is shown in Appendix B. The intersection spectrum 
is defined, as can be seen from Eq. (|3.6|) , as the set of values of E for which Pq{E) lies in 
the range 

± 4tnr.(|) - T,(|)} - {-imti, + fj . (3.16) 

As tb approaches ta this range approaches zero as 

4(t,-t,)t^-^T;(|), (3.17) 
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and each band centered on Ea has width approximately equal to this range divided by 
\PQ{Ea)\, SO that the sum rule for the derivatives of P at the points of the intersection 
spectrum at the critical point ta = tf, > t^i + tab is 

^ 1 1 
For q large and tb = ta > t^i + tab this gives 



<L 1 xtl-Af^ 2^/t2_4t2 



For tb = ta = tab + tab it givcs 



1 ^af) tab 



for all q. For the bicritical case tf) — ta — — 2^^^ it gives 



For tfe > t^i + tafe > there is no intersection spectrum, since the ranges of the constant 
in Eq. (|3.6| ) are nonoverlapping for qk2 = and qk2 = vr. However, we can get an exact 
expression for the intersection over ki of the spectra for fixed ki. This would be the inter- 
section spectrum for the dual problem with ta, t^ interchanged. This spectrum is the set of 
values of E for which Po{E) lies in the range 

± - Kb - Kb] - (-l)^4t%(|) , (3.22) 

and it is shown in Appendix B that the measure of this spectrum is exactly 4(tfe — tai — tab)- 
The same argument that led to Eq. ( p.l8| ) leads without approximation to Eq. ( p.20| ) in the 



case tb = tab + tab > ta, and to Eq. ( p.21| ) in the bicritical case tf, = 2t^i = 2tab > ta- 

These sum rules can be used both to generate rough estimates for the measure of the 
spectrum (the union over k2 of the spectrum for fixed ^2) and to get rigorous lower and 
upper bounds, by repeating the arguments used by Last and Wilkinson0 and Last0. The 
rough estimate of this sum of band widths is obtained by multiplying the range of the 
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constant term in the expression ( p.6|) for P{E) by the appropriate expressions for I] 1/|-P'| 
in Eqs. (CT)-(ra- For h = ta > t^i + tat, this gives 



W ^ 8^tl - Atyq . (3.23) 

For tb = t^i + tab > ta it gives 

W^8\tab-tab\/q, (3.24) 

while for U = 21^1 = 2tab > ta it gives 

W ^ 8tb/q^ . (3.25) 

The argument given by Last0 for the upper bound needs no modifications for the case 
we are considering. The result he obtained is that if the spectrum is defined as the set of 
values of E for which 

-hi< PoiE) < 62 , (3.26) 

where bi, 62 are positive, with Po{E) a polynomial whose zeros E^ are all real and distinct, 
and for which the zeros of the derivative all lie outside the bands, then the sum of the widths 
of the bands W satisfies 

W^<e(6i+62)EW^- (3-27) 

The argument for the lower bound needs some modification because bi 7^ &2j but a straight- 
forward extension of Last's argument gives 

1 



W>- min(6i + v'fof + 46162, &2 + ^Jbl + Ab^b^) J2 T^iT^ ■ (3-28) 

The parameters bi, 62 are given by the differences between the value of the expression in 
Eq. ( p.l6| ) or eq. ( |3.22| ) that defines the intersection spectrum, and the two similar expressions 
that define the band edge. For tb = ta > tab + ^afe we get 

?^i,2 = 8t%(|)±4(t^, + t^,), (3.29) 
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and the second term in this expression becomes neghgible for large enough q. The bounds 
are therefore 



Se^tl - Atyq >W> 2(^5 + l)^tl - Af^/q . (3.30) 

For tb = tab + tab > ta the parameters are 

6i,2 = 4t^, + 4t^,±8t%(|). (3.31) 

For t^i 7^ tab only one of these terms is relevant in the large q limit, so that the bounds 
obtained from Eqs. ( jOq ), (lOTi ), QO^ ) and ( p3TD are 



8e|t,5-U|/g> W^>2(v^+l)| l/g . (3.32) 

For tb = 2tab = 2tab > ta the ta dependent term in Eq. (|3.31|) is also relevant, and 
Eqs. ([Oil , ( F^ , ( lOSi ) and (g) give 



^ > > ^ V^[V^ + ^] ^ (3-33) 

where 

7=|T,(-^)| (3.34) 

This is an example of the importance of the second relevant length scale mentioned in 
connection with Eq. (p.9|). 

For the case ta = tb = 2tab = '^tab this line of argument gives us no useful lower bound, 
since all our lower bounds reduce to zero. For tab + tab > tb > ta we have not succeeded in 
finding an exact expression for the intersection spectrum or some equivalent spectrum. 

IV. MULTIFRACTAL ANALYSIS 

In this Section, we perform numerical scaling analyses for the energy spectrum when 
4> = p/q approaches the quadratic irrational number l/r = — l)/2: the inverse of the 
golden mean. We take p/q to be where F„ is the rath Fibonacci number defined 
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recursively by F„ = + F„_2 and Fq = Fi = 1. Note that g = F„ ~ r" for a large n. In 
order to obtain a spectrum for a rational approximant, the Bloch theorem is applied first to 
Eq. ( |3.1| ), which is periodic with period Then the system becomes effectively finite and 
the spectrum is obtained by a numerical diagonalization. 

To discuss localization of the eigenstates of Eq. ( |3.1| ) we examine its spectrum for fixed 
/c2- When (p = p/q, the spectrum consists of q bands whose widths are denoted by Aj 
{i = l,...,q). Since each band has the same number of states, we assign a probability 
measure 1/q to each band. With increasing q, each band splits into many sub-bands. In 
order to understand the scaling of the spectrum, we introduce a scaling index a by 

l/g-Af. (4.1) 

If states are localized in the limit of g — oo, the band widths A decrease exponentially with 
q, so that a goes to zero. This corresponds to a point spectrum. If states are extended, on 
the other hand, A scales as 1/g, so a is 1. This corresponds to an absolutely continuous 
spectrum. In the critical case, a is expected to take values between and 1. The values 
of a have a distribution on the whole spectrum. This situation corresponds to a singular 
continuous spectrum. 

It is clear from the analysis of the difference between the union spectrum and the in- 
tersection spectrum given in Sec. II10 that the value of ^2 only plays a crucial role for the 
discrete spectrum, when the eigenstates are localized. For the absolutely continuous spec- 
trum there is only an exponentially small dependence on the value of For the critical 
case, since the intersection spectrum of zero measure divides the union spectrum into two 
equal halves, it is clear that to take fixed k2 gives a spectrum whose measure is roughly 
half that of the union spectrum. In the subsequent discussion we will avoid the region of 
localized states, and consider the spectrum found by taking the union over all values of ^2. 

For a systematic analysis of systems with such complex scaling behavior, it is convenient 
to use the multifractal technique developed by Halsey et a/.ii They have introduced the 
spectrum of singularity f{a) defined by 
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n{a) ~ (A)-^(") , (4.2) 

where Q{a)da is a number of bands whose scahng index a hes between a and a + da, and 
(A) is a representative value of A which was not specified clearly in Ref. 23. We first explain 
the multifractal technique as reformulated by Kohmoto.0 

A. Formulation 

First introduce a scaling index e by 

ei = --lnA, . (4.3) 

n 

It is related to a by 

ae = lnr. (4.4) 
We also define an "entropy function" S{e) by 

5(e) = -lnfi'(e), (4.5) 

n 

where fl'{e)de is the number of bands whose scaling index lies between e and e + de, namely 
fi'(e) = Q{a)\da/de\. Here it is important to notice that Aj and ^'{e) depend exponentially 
on n. A band at the nth level splits into many bands at a higher level and may thus yield 
a number of different values of the scaling indices e. However, we expect that the entropy 
function which represents the distribution of e will converge to a smooth limiting form as 
n tends to infinity, and give complete information about the scaling behavior. As in the 
formalism of statistical mechanics, it is convenient to introduce a "partition function" and 
a "free energy" which are defined by 



i=l 



Z„(/5) = EAf (4.6) 

and 



F{P)= hm -InZniP). (4.7) 

n— >oo n 
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Once the free energy is calculated, the entropy function is obtained by a Legendre transfor- 
mation, 

S{e) = F{(3) + (3e, (4.8) 

Thus by changing "temperature" jS one can pick a value of e and then the corresponding 
^(e) is calculated. On the other hand, (3 can be written in terms of e as 

,4.10, 

Usually S{e) is defined on an interval [cj^jj^, emax] and there is no scaling behavior corre- 
sponding to e which is outside the interval and ^(e) = 0. However, F{P) is still defined 
there and from (|4.8|) it is given by F{P) = — emax/5 for /3 > P-^in ^iP) = "^min/^ 
P < /Smax- Thus useful information is only contained in F{P) for the region between 
and /?max where it is not linear. 

Using Eq.( |4.4|) and identifying (A) = exp(— ne) (see (|4.1| )), /(a) is related to the entropy 



function by 



/(«) = (4.11) 



/(a) is defined on an interval [a^[^, amax] where a^un = Inr/emax and amax = ^"^T/^min- 
The maximum value of /(a) gives the Hausdorff dimension of the spectrum. 

B. Numerical results for /(a) 

In all our numerical work we have taken t^i = tab, so we refer just to the value of tab in 
the rest of this Section. Since we are studying the union of the spectrum over all values of 
^2, the system is symmetric with respect to an interchange of ta and tb. Thus we carry out 
numerical calculations only for tb < ta- The system is characterized by two parameters tb/ta 
and tab/ta- 



Before going to the numerical results we notice that the results of Sec. Ill for the Lya- 
punov exponents (in particular Eq. and the discussions below it) determine types 

of spectra in some parts of the parameter space. In region I in Fig. |1] {tab/ta < 1/2), the 
eigenstates are extended in the x direction for irrational (p. Thus the spectrum is absolutely 
continuous and the minimum of a is Oniin — ^ /('^min) — ^- ^i^^ (^a — '^b), 
Eq. (p.l|) is self-dual and the Lyapunov exponents for both x and Px directions are zero. 
Then the spectrum is expected to be singular continuous. The properties of region II have 
not been determined by the analysis of the Lyapunov exponents in Sec. Ill, but it will be 
shown numerically that the spectrum is singular continuous in the whole region II. 

The energy spectra for ta = (self-dual line BD) and n = 10 are shown in Fig. H(a). 
As is well known, the spectrum for tat = (critical point of pure Harper's equation) has 
a self-similar structure. The whole spectrum has three main bands, each main band has 
three subbands, and so on. This structure remains unchanged for < 2tab < ta, but at 
2ta6 = ta it changes. For example, band edges a and b in pure Harper's equation in Fig. |^(a) 
continuously change to a' and b' as tab is increased. They are still band edges for 2tab < ta- 
When 2tab — ^ ta, however, the two points tend to the same point c and are no longer band 
edges. See Fig. 0(b). For 2tab > ta, the topological structure changes further. 

In order to see the global scaling behavior, we have performed numerical calculations of 
f{a), which is defined in the limit of n — oo, by extrapolating the numerical data up to 
n = l9 {q = Fi9 = 6765). 

Figure § is a plot of /(«) for tab = -4. This plot is identical, within the precision of the 
plot, to the plot obtained for the case tab = (pure Harper's equation). As was pointed 
out by Tang and Kohmoto0, /(a) is defined on the interval between a^in — 0.421 and 
«max — 0.547, and takes the maximum value / = 0.5 at a = 0.5. Calculations for a variety 
of other values of tab/ta in the range < tab/ta < 0.5 give results that look identical to 
Fig. 1^, and we conclude that /(a) is universal in this range. At the critical point of pure 
Harper's equation, aj^in i^ the scaling index of the edges of the spectrum (and the edges of 
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each subband), whereas amax is the scahng index of the center of the spectrum (and the 
center of each subband). Even if tab is non-zero, the situation remains the same as long 
as 2tab < ta- More specifically, the scaling index of each band is identical to that of the 
topologically corresponding band of pure Harper's equation. 

At the bicritical point C {2tab = ta = tb) the shape of /(a) suddenly changes. The 
nonzero values of / are found in the range between a-^\^ ~ 0.272 , amax — 0.421, and the 
maximum value of / is about 0.33 at a ~ 0.33. Note that amax at this point is identical to 
'^min 2tafe < ta- At 2tab = ta, howcvcr, the topological structure of the spectrum and the 
scaling indices are different. For example, the scaling index of the bands coming from the 
centers of the subbands of pure Harper's equation is 0.303. It is 0.289 at the edges of bands 
{e.g., d, e and / in Fig. 0(b)). On the other hand, the index remains 0.421 at c in Fig. |^ 
and becomes amax- 

Plots of /(a) on line CD ( 2tab > ta = tb) were obtained for some higher values of tab/ta- 
Although /(a) is not universal on line CD, there is some similarity between the curves we 
found, in that the maximum points are nearly at a = 0.4 and the maximum values are also 
about 0.4. The values of a-^i-^, amax, «0; the position of the maximum of /, and /(ao), the 
maximum value, are shown in Table I. 

We have also calculated the form of / on the line AC {2tah/ta = l,ta > tb). As in the 
previous case on line CD, the curves are somewhat similar to one another in the sense that 
the maximum points are at 0.3 < a < 0.35 and the maximum values are also between 0.3 
and 0.35. The main features of these results are also given in Table I. 

Finally we have calculated /(a) in region II, and one example is given in Table I. In 
this region, the convergence of the extrapolation n oo from the numerical data for finite 
n's is not so good as the previous cases. Thus we cannot obtain reliable results of /(a) to 
claim some universal features. Even though the errors are rather large, however, it appears 
that /(a) is a continuous function defined on a finite range of a. Thus we conclude that the 
spectrum is multifractal and singular continuous in region II. 
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V. NUMERICAL RESULTS FOR BAND WIDTHS 
A. Total band width for Fibonacci sequence 



In the previous Section, we have obtained the following: (1) for ta = and 2tab < ta (line 
BC in Fig. |l]), the scaling behavior of the spectrum is universal, that is, /(a) is completely 
identical to that of the tab = case ( pure Harper's equation); (2) at the bicritical point of 
ta = tb and 2tab = ta (point C in Fig. |1|), the scaling behavior suddenly changes; (3) when 
ta = tb and 2tab > ta (line CD in Fig. |I]), the scaling behavior is clearly different from those 
of (1) and (2). Although it is not completely universal, the changes are small for increasing 
tab/ta'i (4) when ta > tb and 2tab = ta (line AC in Fig. |l]), the scaling behavior is similar 
but not identical to that of point C {ta = tb and 2tab = ta)] (5) when 2tab > ta and ta > tb 
(region II in Fig. |l]), it is difficult to estimate f{a) but it is certain that the spectrum is also 
multifractal and singular continuous. 

To make the above statements more concrete, we investigate scaling of the total band 
widths. Recall that in pure Harper's equation in the critical case (ta = tb,tab = 0), the total 
band width W scales as W ^ 1/q for large q (see Eq. (|1.5|) ). We know from the analytic 
results of Sec. HI that this scaling also holds all along the line BC (ta = tb > 2tab) except 
at C. The numerical results show that for the Fibonacci sequence the result 



holds accurately, in agreement with Eq. ( p.l4| ). 

On the bicritical line AC {2tab = ta > tb), which separates the region of ta dominant 
from the region of tab dominant, we know from the results of Sec. HI that the sum of the 
band widths should scale as q~^. Figure ^ shows the results for q'^W plotted against n 
(~ Ing/lnr) for different values of tb/ta on this line. In this subsection all the plots of W 
are in the units such that ta = 1. It is seen that the scaling index 6 for the global behavior 
of W which is defined by 




(5.1) 



W ~ (l/q)' , 



(5.2) 
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is 2. For ta = h = 2tab, the quantity WF^/ta tends rapidly to 6.4911. In addition to this 
power-law decrease, an oscillatory behavior is observed for 7^ ta- The oscillation has a 
period 4 against n in the case 2tb = ta, which is probably related to the period 8 of the 
Fibonacci sequence modulo 3. 

In the region 2tab > ta > tf, we find two different types of behavior, but in all cases the 
analysis is complicated by oscillatory terms superposed on a general power law dependence 
on q. Plots of \nW against lnF„ (which is proportional to n) appear to lie close to a line of 
slope —1.25 for the case 2tab > ta = tb (the line CD in Fig. |lD, whereas they cluster around a 
line of slope —1.56 for the case 2tab > ta > tb (the interior of the region II in Fig. p. Figure]^ 
shows plots of \n{F^-'^^W) against n for various examples of 2tab > ta = tb- There is a period 
4 oscillation in the case tab = ta = tb, similar to the oscillation of period 4 that shows up 
in Fig. ^ for the case tab = taf^ = tb- There is also an oscillation of period 6 for the case 
tab = ta/ V2. In other cases the oscillation around the general horizontal trend is irregular, 
and show no signs of diminishing as n increases. An accurate estimate of the exponent 6 
can be made when the period is short, but we can only make a rough estimate when there 
is no period within the range of n we can use. There appears to be a slight difference, of 
order .01, between the values we get for S with tab = ta = tb and with tab = ta/ \f2 = tf,/^/2. 

Figure ^ shows plots of \n{F^-^^W) against n for various examples of 2tab > ta > tb. 
Again, there are large oscillations about the general horizontal trend of the plots, and a 
simple period, 6 in this case, can be seen clearly for tab = ta/V2, t^ = 0. For reasons which 
are discussed later in this paper, we think that this period is related to the period 12 of the 
Fibonacci sequence modulo 4. The value of 5 = 1.56 is not so easy to estimate for these 
examples, but it is clearly intermediate between the value we found on the line CD of Fig. |l|, 
and the value of 2 which is known to be correct for the bicritical line AC. 

The results in this subsection are summarized as follows. The scaling index 6 for the 
total band width is 1 on the critical line BC in Fig. |^. Not only at point C but also on line 
AC, 6 is 2, i.e., the bicritical behavior is found. On line CD, we cannot find a significant 
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dependence of S on tab/ta, and it is about 1.25. Also in region II, the situation is similar 
and the value of 6 is about 1.56. 

In the limit of tab/ta —>■ oo, the problem reduces to the case only with nearest neighbor 
couplings. Thus the 6 must be 1 in this limit. Although we investigated 6 for quite large 
tab/ta (up to tab/ta = 10) in rcgiou II and on the line CD, we could not find a tendency that 
6 decreases and approaches 1. 

B. Scaling for (p = l/q 

For rational approximants to the golden mean the total band width is spread over a large 
number of bands in different energy ranges. For large denominator rational approximants 
to a small denominator rational, say po/qo, the total band width is concentrated in narrow 
ranges about the go values of the energy where there is a logarithmic singularity in the 
density of states for the case = po/o'o-^'^ For sequences such as0=l/gor2/g the limit 
of the sequence gives po = 0, qo = 1, and all the width comes from one singular energy at 
or near the center of the band. To a considerable extent the results for this case can be 
understood in terms of the WKB analysis presented in Sec. VI, and many, but not all, of 
the results appear to generalize to more complicated sequences of fractional values of (p. 

We have made some numerical checks of the scaling relations ( |1.5| ) and ( 3.13| ) for the case 



tb > ta > '^tab, and of the scaling relation (|3.15|) for tb > 2tab > ta, but these have not been 
extensive, and we have not found particularly interesting features. We have concentrated on 
three critical or bicritical cases. For tb = ta > 2tab the bounds ( ^.3(J| ) show that the width 
must scale as 1/g, and finite size scaling theory suggests that it should have the limiting 
form given by Eq. ( |3.14| ). For the bicritical case tb = 2tab > ta the bounds Eq. ( ^.33| ) show 



that the width must scale as 1/g^ (except possibly in the special case ta = tb), but Eq. (FTSD 



does not tell us much more, since we expect the additional length given by Eq. (|3.9|) to be 
involved. Finally, there is the critical case 2tab = 2tab > tb >ta, which one might expect to 
be analogous to the other critical case, since the dominant terms are essentially the same. 
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but rotated through an angle 7r/4 in the Px, x plane, with half the size of unit cell. However, 
we know no rigorous bounds in this case, have derived no useful finite size scaling relations, 
and expect the influence of the additional lengths given by Eq. ( p.9| ) to be important. 

The case ta = h > 2tab is critical and the measure W scales like \Jtl — 4t^^/ q as predicted 
by the finite size scaling theory of Eq. ( p.l4| ). We have done a numerical check for numerator 
2 and q odd, where we have evaluated qW/ \Jt'i — 4t^^ for tab=^, 0.2, 0.4 and found that they 
very rapidly converge to a common value. The difference between values of qW/ \Jtl — 4t^j 
for different values of tab is less than one part in 10^ for q greater than 41. This implies 
that the energy scale is reduced by a factor y^l — 4t^^/ 1"^ and that there are no logarithmic 
corrections in the case of p=2 as we increase tab- 

For numerator equal to unity, the scaling limit remains the same, but there are large 
corrections to scaling, which are shown in Fig. |^. These show an interesting cusp-like 
oscillation of qW/ y^l — At^^/tl for non-zero tab- The periodicity of this oscillation can be 
related to the integral of the classical momentum over the length of the system, as we discuss 
in Sec. VI, while the magnitude of the oscillation is bounded by two curves which are followed 
by odd and even values of q with ^^{,=0.1 

The 2tab > max(ta, tb) region is analogous to the ta = tb dominant region just considered 
in that, without ta,tb, the problem reduces to that of the Harper's equation with twice as 
much flux per unit cell, and that ta, tb may be regarded as perturbations to the critical Harper 
problem. It is probably not important that the perturbations generally break the square 
symmetry of the system, since, as we discussed in Sec. II, the symmetry under —px 
remains, as well as rotation of phase space by vr. In the cases we have examined with = 1/g 
the measure scales like 1/q for ma.x(ta,tb) < 2tab- As was mentioned in Sec. HI, qW does 
not seem to tend to a limit, but oscillates in the tab dominant regime. 

To study the periodic nature of the measure, it is convenient to choose the set of param- 
eters ta/2tab = cos(7rpi/gi), tb/2tab = cos{7ip2/q2), and characterize the system with a set of 
fractions {pi/qi,P2/(l2)- Several graphs of qW as a function of q are shown in Fig. p. For 
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Pi = P2 = 1, we have found strong peaks at multiples of qi x q2 (primary peaks) and much 
less strong peaks at multiples of qi (secondary peaks) if qi is a small integer such as 2 or 3. 
For qi and q2 both large, close, and relatively prime, for example, 5 and 6 or 5 and 7, the 
secondary peaks do not seem to occur at definite multiples of either number. The primary 
peaks always occur at multiples of qi x q2. If qi = q2 or they share a common factor, then 
the secondary peaks occur at integer multiples of their lowest common multiple, and the 
values of qW at these points nearly match those of the primary peaks. 

In the bicritical case with 2tab = tb > ta we know from the inequalities (|3.32| ) that the 



measure of the spectrum must be proportional to g ^, and there is one length scale remaining 
which is given by Eq. ( |3.9|) , so we might expect a scaling form 

q'^W = ghc{qaTccos{ta/2tab)) ■ (5.3) 

Figure]^ shows gbc{x), with x = q aYccos(ta/2tab) , for ta/2tab=l/S and 0.99. In the limit ta = 
'^tab, gbc{x) becomes a slowly increasing function of x, and we have not determined its limiting 
value; this slow convergence may be a special feature of particular fractional forms of (p, as we 
found no sign of it for the Fibonacci sequence. Figure ^(b) transparently displays a cusp-hke 
form of the scaling function, whose variation lies well within the bounds given by Eq. ( ^.33[ ). 



Earlier in this subsection, such cusps were reported when we considered the ta = tb dominant 
region. It turns out that a function of the form gbc{x) = A — B ln(l + | sin(7r(x — 6)\), where 
A, B and 5 are constants chosen to fit, describes the actual scaling function rather well. The 
motivation for this form of gbc{x) came from the idea that x must be a dimensionless quantity 
and that such a quantity can be obtained by multiplying q on both sides of Eq. (|3.10|) before 
one takes the limit g — oo. Since taj, = tab, this gives us x = qX^ = ln2 — ln(H- \Tg(ta/2tab)\) 
which seems to contain essential features of gbc{x). Based on this hypothesis, the cusp-like 
behavior can be understood as the refiection of the importance of the absolute value of the 
Chebyshev polynomial. 
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C. Scaling for (f> = p/{p^ ± 1) 



The sequences = 'pjijP' ± 1) are intermediate between ip = 1/q and = -F„/F„_i. 
The continued fraction expansion has only two terms in it, rather than the n/2 terms for 
the Fibonacci sequence, and it represents a slow modulation of the diagonal term, and so 
could be treated by WKB methods (although we have not succeeded in carrying out such 
an analysis in this case). It is, however, more complicated than the = 1/g case in that 
the significant contributions to the band width come not only from bands in the immediate 
vicinity of the singular energy —tatb/tab, but also from the centers of neighboring clusters of 
p subbands. 

We have studied the bicritical case tb = 2tab=l and ta = cos{piTc / qi) where qi was kept 
small (< 7). With this choice of ta we expect behavior periodic in q in the limiting values of 
q'^W coming from the constant term Eq. (|3.6| ), but there may also be some dependence on 
the numerator p coming from other terms in the characteristic equation. We have confined 
ourselves to numerators not exceeding 20, and denominators up to 401. As in the previous 
case of a simple fraction of 0, the scaling function showed periodicity (up to corrections to 
scaling) with periods equal to qi for both ± 1. The pattern of graphs are quite different in 
two cases and values of maxima and minima as well as where the maxima and minima occur 
do not agree in general. The convergence to the limit was much slower here than it was with 
= 1/q. This might have to do with the fact that, for the fraction = p/{p'^ ± 1) ~ l/j9 
with p < 20, the corrections are not yet completely negligible. For pi/qi = 1/4, we should 
expect a period of 2 because ± 1 (mod 4) alternates between 2(0) and 1(3) for odd and 
even p but instead we observe a period of 4. Apparently the scaling function here is more 
complicated than what was the case if the relevant variable was simply the ratio of two 
length scales of the system. This shows conclusively that there is periodic dependence on 
the value of the numerator as well as on the value of the denominator. 

In the tab dominant critical region we have found instances of an anomalous (non-integer) 



scaling exponent. Figure [1^ shows a plot of \og{qW) against log q for the case 2tab = 1, = 0, 
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tb = cos(7r/4), = + 1). For numerators equal to 2 mod 4 we see the points lying 

on a line with a negative slope, whereas linear scaling should give no slope at all, as we 
see for the other values of p. The slope both of this plot, and of the very similar plot for 
(f) = p/Ip"^ — 1), is found to be within 0.3% of -1/2, which implies that, for the sequence 
p/q = 2{2n + l)/[4(2n + 1)^ ± 1], the scaling is like l/g^/^ rather than l/q. For = 0,4 = 
cos(7r/3), the plot alternates between points with exponent close to 1, if p = 0, ±1 (mod 6), 
and those with exponent close to 1.5 if p = 2, 3, 4 (mod 6) for both types of denominator. 
It is much harder to extract the critical exponents here because we have to increase p by 6 
instead of 4 to arrive at points lying on the same line. For {] < ta = tb < 2tab, the evidence for 
non-integer exponent is far less obvious and we are not sure if the exponent is significantly 
different from 1. 

The noninteger exponents found in our studies of the Fibonacci sequence are likely to be 
related to these results. 

VI. WKB THEORY 

The saddle point value of the energy contour Eg corresponds to quantum mechanical 
states that can thread the system without attenuation and there exists an interesting relation 
between the integral of the classical momentum p^ given by Eqs. (|2.1| )-( p^ between turning 
points with periodicities in the scaling functions. 

In Ref. ^ it was shown that for = 1/g the sum of the band widths could be related to 
the Green function at the turning points. For ta = tb > tai + tab, energy close to —2(tg^b + tab) 
and X = (j)/2 + C,, where ^ is small, the continuum approximation for Eqs. ( |1.6| ) and ( p.l|) 
takes the form 

H + 2{t,b + tab) ~ 

A7T^(f)\ta + t,b + tab)-^ + (4 " t,b ' tab)i^ + 27r#(t,-, - tab){i^ + ^0 ■ (6.1) 
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A very similar expression is obtained near the other turning point with x close to zero and 
Px close to TT. This can be diagonalized by a canonical transformation, and the energy scale 
it yields is 



47r0^t2 _ ^t^.tab (6.2) 

in both cases. The earlier arg umentsiiffl applied to this case give the scaling result quoted 
in Eq. (|3l4l ). 



The corrections to this scaling form can be calculated by using the WKB approximation 
to get a more precise approximation to the band width^, which involves the connection 
between the turning points as well as the behavior at the turning points. We have not 
worked this case out in detail, but we know that for 2tgj, = 2tab the phase change around an 
orbit close to the critical orbit is 

, - I , I 2tab + taCOSX\ 

p^ax — / dxarccos . (,6.oJ 



27r0 / TT(j) Jo \ ta + 2tab cos X, 

Differentiation of the right side of this equation with respect to tab gives an integral that 
can be evaluated explicitly, and reintegration of this result gives 

^ i pjx ^ gvr - ^ / ^^d9 . (6.4) 

2n J n Jo sinh 9 

The first term in this expression gives rise to the correction to scaling that alternates with 
the parity of q even for tab = 0, while the second term gives a period in q that goes like 
7!''^ta/4:tab for Small tab- FoT tab/ta = 0.1 it givcs 24.5 as the period, which is in good agreement 
with the numerical results shown in Fig. 0. In the limit 2tab — ^ ta the second term is vrg and 
cancels the periodicity due to the first term. 

For 2tab = 2tab dominant the problem is somewhat different. The classical contours 
through the saddle points are given by 

-^{tb + 2tabCOSp^){ta + 2tabCOSx) = . (6.5) 
tab 

The quadratic approximation to the Hamiltonian near one of the saddle points is given by 
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H H ^ ± -I, 



+ (6.6) 



and diagonalization of this by a canonical transformation gives an energy scale 



This quantity gives a good account of the relative sizes of the energy scales of the band 
widths shown in Fig. |. However, as we showed in the study of the characteristic equation 
in Sec. Ill, the lengths given in Eq. (|3.9| ) are certainly relevant, and there must be terms 
periodic or nearly periodic in q. These should come from the rectangular contours given in 
Eq. ( |6.5|) , whose areas are 

4 arccos(±— ^) arccos(± — ^) , (6.8) 



and it is the ratio of these areas to the quantum of action An'^cj) given by Eq. (|1.2|) that 
determines this periodicity. 

Qualitatively this accounts for the periods in g of 9 and 21 which show up in Figs, j^a) 
and (b), since the smallest areas given by Eq. ( |6.8| ) are 471^/9 and 477^/21 in the two cases, 
and the larger rectangles are multiples of these. To understand these results in more detail 
we need to make a more careful study of the way the phases affect the dynamics. 

The bicritical case th = 2tab is simpler, since one of the arccosines in Eq. (|6.8|) is equal 
to TT. For ta = tab the two rectangles have areas 471^/3 and Svr^/S, so the period 3 in g which 
can be seen in Fig. H(a) should be expected, while for ta = A95tab, vr/ arccos(ta/2tafe) is equal 
to 22.2, which agrees well with the period shown in Fig. ^(b). 

For a large denominator fraction p/q approximating a small denominator rational po/qo, 
the commutator remains finite and the WKB approach does not directly apply. 

Wilkinsonii has shown that at ta = tb,tab = tai = 0, each one of go clusters of bands 
are described by an effective Hamiltonian obtained from quantizing the inverse of the char- 
acteristic polynomial (|3.6|) for the nth band. 



En{h, k2) = Po:n(2C cos(goA;i) + 2^ cos(goA;2)) . (6.9) 
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The new flux for each cluster depends on the Chern number for that particular band, but 
is in general of the order of the difference p/ q — Po/ qoi therefore small. A more direct waylll 
is to regard the problem as having cf) = Po/qo with the wavevector /c2 slowly modulated with 
a period qsgo/hsPo -PsQol- The /c2 in (|3^ ) becomes /c2,n = k2 + 27m{{psqo -Poqs)/qoqs) and 
ki is also site-dependent. This approach generalizes to tab 7^ with the result that singular 
energies are the go roots of the equation 

Po{E) - A{-irC = (6.10) 

for ta = tb dominant and of 

Po{E) - A{-iyHZT,,{^)T,,{^) = (6.11) 

for tab = tab dominant regime. We have done some numerical check for po = l,qQ = 2 



and found that the singular energies come at ±2ta and ±2^tl + tl — t^tf respectively in 
agreement with predictions given by Eqs.( |6.lOD and ( |S.ll ). 



VII. DISCUSSION 

Our studies have shown that the characteristic critical properties of Harper's equation 
persist when the symmetry is broken by terms which couple sites to their next nearest 
neighbors, provided that the reflection symmetry in the diagonals is preserved, and provided 
that the next nearest neighbor terms satisfy the inequality tab + tab < ta = tb- In this 
region the multifractal analysis gives a universal result, strict bounds for the width W of 
the spectrum show that it must go to zero with the reciprocal of the denominator q, and 
numerical results and semiclassical analysis show that qW has a limit which is rescaled by 
the next nearest neighbor coupling, but which is independent of the numerator p. It is only 
in the corrections to scaling that any oscillatory behavior shows up. 

For the bicritical line tab = tab = iriax(ta, ^6)/2, which divides the region dominated by 
the nearest neighbor terms from that dominated by the next nearest neighbor terms, the 
multifractal behavior is quite different, and the width W of the spectrum is known to be 
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proportional to as we know from the rigorous bounds,ll2lO as well as from numerical 

work. In this case an oscillatory behavior superposed on the dependence shows up, 
whose periodicity is at least partially understood from the WKB analysis of Sec. VI. 

The region for which we have least understanding is the region dominated by the next 
nearest neighbor terms, with tab = ial > niax(ta, 4)/2. Although the case ta = Q = h is 
equivalent to the original Harper equation, with the axes turned by 7r/4 and the unit cell 
doubled in area, for any nonzero values of the nearest neighbor coupling their strengths 
taitb remain relevant, as can be seen clearly from Eq. |3.6| . The spectrum has a multifractal 
form, but the multifractal analysis gives different ranges of the exponent for different values 
of the parameters of the Hamiltonian. There seems to be a marked difference between the 
behavior for ta = h and for 7^ tf,. When the width W of the spectrum is studied the 
results are quite different for the ratio of Fibonacci numbers and for 1/g or similar sequences 
of fractions with fixed numerator. For = 1/g we find W going to zero like with an 
oscillatory coefficient, while for the Fibonacci sequence the oscillatory behavior is similar, 
but the dependence on the denominator is of the form where 5 seems to be about 1.25 
for ta = tfe, and 1.56 for for ta 7^ tb- 

We believe a partial understanding of this behavior can be obtained from our results for 
of the form l/(gi + 1/ 52) which were discussed in Sec. [V Q . For large values of the numbers 
in the continued fraction expansion of we get relatively simple scaling behavior which can 
be studied by using WKB theory, although we have not carried out the analysis in detail 
yet. Because of the periodic or nearly periodic behavior as the qj are varied, each stage of 
the scaling may carry one arbitrarily close to the bicritical boundary of the critical region. 



and we found examples, one of which is shown in Fig. |T0|, where for certain values of qi the 
dependence of W on q2 at the next stage of scaling was the l/g| typical of the bicritical 
boundary. Approximants of the golden mean which are the ratios of two Fibonacci numbers 
have continued fraction expansions in which the terms are particularly small, so one is very 
far from the simple scaling behavior expected for large qj, so the scaling behavior at each 
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stage may be intermediate between the critical and the bicritical behavior. 
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APPENDIX A: EVALUATION OF PRODUCT OF OFF-DIAGONAL TERMS 

Evaluation of the product in Eq. ( |3.3| ) is required to obtain Eq. ( |3.6| ). The expression 



P{k,) = n (^a + hie-^^""^--^-''' + ta.e^("+^)+^'=^) (Al) 



n=0 

must be periodic in ^2 with period 27r/g, since the addition of 27r/g to ^2 yields a permutation 
of the same factors in the product. This tells us that the only terms that survive in the 
product are those with q factors of exp(— 2/^2), those with q factors of exp(z/c2), and those 
with equal numbers of factors of exp(— zA;2) and exp(ifc2)- The first two cases give 

^q_^-TTipq-iqk2 _|_ ^^nipq+iqk2 (A2) 

while the k2 independent terms can be written as 

9-1 



Q=Yl{ta + ^ ^g— _ (_i)P92t9 cos(gA;2) , (A3) 

n=0 



where = taitab] here we have subtracted off the k2 dependent terms using Eq. ([A2|) . The 
expression obtained from Eq. ( [A3[ ) by setting ^2 = 0, 



Q/t" + {-Vfn = n (ta/t + e- — ("+^) + e— ("+^)) , (A4) 



{ta/t + e 

n=0 

is a polynomial in ta/t whose zeros are given by 



y = -2cos[^(n + i)]. (A5) 



These are the zeros of the equation 
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cos(garccos|) ^ T,(|) = {-If-" = -{-ir , (A6) 

where Tg is the Chebyshev polynomial of order q. Since the coefficient of x'^ in Tg{x) is 2"^"^, 
this gives 

Q/ti + {-iyi2 = 2Tg{tj2t) + (-1)P«2 . (A7) 
Combination of this with Eq. in Eq. ([Al| ) gives 

n=0 



which is the result used to derive Eq. (|3.ti|) 



APPENDIX B: EXACT EXPRESSIONS FOR THE INTERSECTION SPECTRUM 

Avron, Mouche and Simoni have shown that for the case t^^ = Q = tab the intersection 
spectrum has measure A\tb — ta\. By exploiting the techniques used in earlier wor we can 
generalize this result to the case h > ta > tab + tab- This argument depends in its details 
on the parities of p and q, so we will give the argument explicitly for the case of p, q odd. 
For qk2 = vr and qki zero or vr we exploit the symmetry of Eq. ( |3.1|) about the points n = 
and n = q/2, which allows the problem to be reduced to eigenvalue problems for tridiagonal 
matrices of order (g± l)/2. The eigenvalues E^'^ and which correspond to eigenstates 
of Eqs. (|3.1|) and ( p.2|) that are even about both these symmetry points or odd about both 
of them give solutions with qki = 0, while E^~ and -E""*" give solutions with qki = vr. We 
number the eigenvalues from highest to lowest, starting with zero for the solutions Eq^ and 
Eq^ which are even about n = 0, and starting with unity for those that are odd about this 
point. The highest value of m is (g — l)/2 in all four cases. The matrices corresponding 
to odd and even boundary conditions about the point n = differ only in whether the 01 
element is zero or not, and the trace formula gives 
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24 - E+- + J2{E-- - E+-) = , (Bl) 

m 

where each term in the sum is positive. The matrices corresponding to odd and even 
boundary conditions about n = q/2 differ in having an extra term ±{ta — t^j, — tab) in the 
lowest diagonal element, so the trace formula gives 

~ -^m") = 2(ta - ia6 - ^af>) , (B2) 

m 

with each term in the sum positive. Addition of these two equations gives 

- E+-) = Et - 2ih -ta + + U) . (B3) 

m 

The terms in the sum are all positive, and are equal to the band gaps with qk2 = 0, 
qki = vr in the intersection spectrum, while the eigenvalue Eq~ is the highest band edge 
of the intersection spectrum. A similar argument can be constructed for the solutions with 
qk2 = IT, qki = 0, where we number the eigenvalues S^'^, in increasing order. This gives 

++ - = -S^+ - 2{U -ta- - . (B4) 

m 

Again, all the terms in the sum are positive and give band gaps of the intersection spectrum, 
while £q'^ is the lowest band edge of this spectrum. Addition of Eqs. (p3|) and ( p4D gives 
the measure of the intersection spectrum as 



For the case U > t^i + tab > ta there is no change in the argument that leads to Eq. (|B4D , 
but Eqs. (pT|)-(p3|) must be replaced by 



2tb - E++ + - = , (B6) 

m 

- = 2{tal + tab - ta) , (B7) 



and 
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- -^m"^) = - 2(tb + ta- t^l - tab) , (B8) 

m 

with all terms in the sums positive. Addition of this to Eq. ( [B4| ) gives the measure of the 
spectrum for qki = as 

m m 

and this gives the generalization of the result for the intersection spectrum for the case 

tb > + tab > ta- 



35 



REFERENCES 



^M. Ya. Azbel, Zh. Eksp. Teor. Fiz. B46, 929 (1964) [Sov. Phys. JETP 19, 634 (1964)]. 

^D. R. Hofstadter, Phys. Rev. B14, 2239 (1976). 

^S. Aubry and G. Andre, Ann. Isr. Phys. Soc. 3, 133 (1980). 

^D. J. Thouless, Phys. Rev. B28, 4272 (1983). 

^See B. Simon and B. Souillard in Random Media ( Springer- Verlag, 1987), IMA Vol. 7, 
and the references therein for a rigorous discussion of the relation between the Lyapunov 
exponent and the localization of eigenf unctions. 

^D. J. Thouless, Commun. Math. Phys. 127, 187 (1990). 

^ J. Avron, P. H. M. v. Mouche, and B. Simon, Commun. Math. Phys. 132, 103 (1990). 

«D. J. Thouless and Y. Tan, J. Phys. A24, 4055 (1991). 

^G. I. Watson, J. Phys. A24, 4999 (1991). 

°Y. Last and M. Wilkinson, J. Phys. A25, 6123 (1992). 

^ Y. Last, Commun. Math. Phys. to be pubhshed. 

^D. J. Thouless and Y. Tan, Physica A177, 567 (1991). 

^M. Kohmoto, Phys. Rev. Lett. 51, 1198 (1983). 

^C. Tang and M. Kohmoto, Phys. Rev. B34, 2041 (1986). 

^S. C. Bell and R. B. Stinchcombe, J. Phys. A22, 717(1989). 

^I. M. Suslov, Zh. Eksp. Teor. Fiz. 83, 1079 (1982) [Sov. Phys.-JETP 56, 612 (1982)]. 

^C. M. Soukouklis and E. N. Economou, Phys. Rev. Lett. 48, 1043 (1982). 

^H.Hiramoto and M. Kohmoto, Phys. Rev. Lett. 62, 2714 (1989). 

^H. Hiramoto and M. Kohmoto, Phys. Rev.B40, 8225 (1989). 

'°B. Helffer and J. Sjostrand, Mem. Soc. Math. France, 40, 139 (1989). 

'^F. H. Claro and G. H. Wannier, Phys. Rev. B19, 6068 (1979). 

'2 Y. Hatsugai and M. Kohmoto, Phys. Rev. B42, 8282 (1990). 

T. C. Halsey, M. H. Jensen, L. P. Kadanoff, L. Procaccia and B. I. Shraiman, Phys. Rev. 

ASS, 1141 (1986). 



36 



M. Kohmoto, Phys. Rev. A37, 1345 (1988). 
M. Wilkinson, J. Phys. A20, 4337 (1987). 



37 



FIGURES 



FIG. 1. Phase diagram of the different regions of the parameters of the Hamiltonian, for 
tab = ^ab- region I states are pure extended, and in region III the states are purely locahzed 
states. On the hne BC the Hamiltonian has both diagonals as reflection lines, and the behavior 
is critical, very similar to that at the Harper point B. Bicritical behavior is found not only at the 
point C, but along the lines AC and CE. 

FIG. 2. (a) Spectra for various points on line BD in Fig. [l| {ta = tb)- (b) Enlarged version of (a). 

FIG. 3. (a) Plot of f{a) for tab = 0.4ia, ta = tb- Other points on the line BC of Fig. | give 
identical results, (b) Plot of /(a) for the bicritical point C where tab = 0.5ta, ta = tb- 

FIG. 4. Plot of F^W versus n for tab/ta = 0.5 and various values of tb/ta- A period 4 
oscillation as a function of n for the case tb/ta = 0.5 can be seen clearly superposed on the l/F^ 
dependence of W. 

FIG. 5. Plots of F^''^^W on a logarithmic scale versus n for ta = tb and various values of tab/ta- 
Period 6 can be seen for tab/ta = l/\/2, and period 4 for tab/ta = 1. 

FIG. 6. Plots of F^'^^W on a logarithmic scale versus n for various values of tab/ta and 
ta/tb 7^ 1- For tab/ta = l/\/2 a period of 6 can be clearly seen. 



FIG. 7. Plots of qW/^l — 4i^j as a function of q for tab=0, 0.1 with ta = tb=l. The squares 
representing tafe=0 form envelopes (upper and lower bounds) for other values of tab < 0.5. 



FIG. 8. Plots of qW/ y (1 — t1){l — t^) as a function of q in the critical regime 2tab = 1 > ta,tb 
for several values of {ta,tb) = (cos(7rpi/(7i), cos(7rp2/92))- Each graph is characterized by a set of 
fractions, ipi/qi,P2/q2)- 

FIG. 9. Plots of q'^W as a function of q at the bicritical point 2tab=^=ta for tb=0.5 and 0.99 
respectively. 

FIG. 10. Plot of qW on a logarithmic scale against Inq for (j) = p/q = p/{p'^ + 1). Points for 
p = 2 mod 4 are clearly seen to lie on a line with a slope -0.50. 
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TABLES 



TABLE I. Table of the extremal values of a and the position of the maximum in the multifractal 
analysis for various values of the parameters of the Hamiltonian. 





h/ta 






ao 


/(«o) 


0,.2,4 


1.0 


.421 


.547 


.50 


.50 


.5 


1.0 


.272 


.421 


.33 


.32 


.5 


.75 


.282 


.381 


.35 


.34 


.5 


.5 


.281 


.366 


.34 


.34 


.5 


.25 


.281 


.370 


.33 


.32 


1.0 


1.0 


.300 


.650 


.43 


.41 


2.0 


1.0 


.305 


.640 


.42 


.41 


3.0 


1.0 


.313 


.628 


.44 


.42 


.6 


.5 


.29 


.41 


.36 


.35 
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